# This file re-runs HPBM results
# and our replication of HPBM at agency and county level
# and stores all of the results
library(plm) #for panel data
library(lfe) #for felm function
library(rlang)
library(tidyverse)
library(broom)
library(foreign)
library(purrr)
library(texreg)
library(foreign)
source("lib/model-functions-HPBM.R") #these have all the different models we need to run
options(scipen=0, digits=7) #to standardize number of digits when exporting data. all data here is run using these options.

controlled_subset <- FALSE
year_subset <- FALSE
city_outliers_eliminated <- TRUE
crime_problem_obs_dropped <- FALSE
scale_instruments <- TRUE

# load data
if(controlled_subset){
  load('data/city-aggregate-controlled.RData')
  load('data/county-aggregate-controlled.RData')
  load(file='data/harris.RData') #Harris data for exact replication
  harris <- df
  rm(df)
} else {
  load('data/city-aggregate.RData')
  load('data/county-aggregate.RData')
  load(file='data/harris.RData') #Harris data for exact replication
  harris <- df
  rm(df)  
}

# restrict data
if(year_subset){
  agency <- agency[agency$year >=2010 & agency$year<=2012,]
  county <- county[county$year >=2010 & county$year <=2012,]
  harris <- harris[harris$year>=2010 & harris$year<=2012,]
}

if(city_outliers_eliminated){
  agency <- agency %>% filter(allcrime_rate < 1000000)
}

# scale instruments
if(scale_instruments){
  agency$zHUdI1_lag <- scale(agency$zHUdI1_lag)
  agency$zHUd1_lag <- scale(agency$zHUd1_lag)
  agency$zHUdI6_lag <- scale(agency$zHUdI6_lag)
  agency$zHUd6_lag <- scale(agency$zHUd6_lag)
  agency$zHUdIland_lag <- scale(agency$zHUdIland_lag)
  agency$zHUdland_lag <- scale(agency$zHUdland_lag)
  agency$zHUdIhidta_lag <- scale(agency$zHUdIhidta_lag)
  agency$zHUdhidta_lag <- scale(agency$zHUdhitda_lag)
  
  county$zHUdI1_lag <- scale(county$zHUdI1_lag)
  county$zHUd1_lag <- scale(county$zHUd1_lag)
  county$zHUdI6_lag <- scale(county$zHUdI6_lag)
  county$zHUd6_lag <- scale(county$zHUd6_lag)
  county$zHUdIland_lag <- scale(county$zHUdIland_lag)
  county$zHUdland_lag <- scale(county$zHUdland_lag)
  county$zHUdIhidta_lag <- scale(county$zHUdIhidta_lag)
  county$zHUdhidta_lag <- scale(county$zHUdhitda_lag)
}

#drop problem observations for UCR data
if(crime_problem_obs_dropped){
  county <- county %>% filter(lessthan100coverage==0) #this drops all county obs that have less than 100% coverage from ORIs in that county-year
} 

# create panel data
agency <- pdata.frame(agency,index=c("g","year"))
county <- pdata.frame(county,index=c("g","year"))
harris <- pdata.frame(harris, index = c("fips", "year")) 

#------------------------------------------------------------#
## Harris tables for supplementary materials ----
#------------------------------------------------------------#

## Harris Replication using our Agency level data -----
homicide.items.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'items')
homicide.value.agency.Harris <- run.IV.agency.Harris(Murder_rate, agency, 'value')
robbery.items.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'items')
robbery.value.agency.Harris <- run.IV.agency.Harris(Robbery_rate, agency, 'value')
assault.items.agency.Harris <- run.IV.agency.Harris(Assault_rate, agency, 'items')
assault.value.agency.Harris <- run.IV.agency.Harris(Assault_rate, agency, 'value')
vehicle.items.agency.Harris <- run.IV.agency.Harris(MV_theft_rate, agency, 'items')
vehicle.value.agency.Harris <- run.IV.agency.Harris(MV_theft_rate, agency, 'value')


# reduced form at agency level
homicide.red.agency.items.harris <- reduced.agency.Harris(Murder_rate, agency, 'items')
homicide.red.agency.value.harris <- reduced.agency.Harris(Murder_rate, agency, 'value')
robbery.red.agency.items.harris <- reduced.agency.Harris(Robbery_rate, agency, 'items')
robbery.red.agency.value.harris <- reduced.agency.Harris(Robbery_rate, agency, 'value')
assault.red.agency.items.harris <- reduced.agency.Harris(Assault_rate, agency, 'items')
assault.red.agency.value.harris <- reduced.agency.Harris(Assault_rate, agency, 'value')
vehicle.red.agency.items.harris <- reduced.agency.Harris(MV_theft_rate, agency, 'items')
vehicle.red.agency.value.harris <- reduced.agency.Harris(MV_theft_rate, agency, 'value')

## Harris County Replication -------
homicide.county.items.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'items')
homicide.county.value.harris <- run.IV.county.Harris(MURDER_countysumrate, county, 'value')
robbery.county.items.harris <- run.IV.county.Harris(ROBBERY_countysumrate, county, 'items')
robbery.county.value.harris <- run.IV.county.Harris(ROBBERY_countysumrate, county, 'value')
assault.county.items.harris <- run.IV.county.Harris(AGASSLT_countysumrate, county, 'items')
assault.county.value.harris <- run.IV.county.Harris(AGASSLT_countysumrate, county, 'value')
vehicle.county.items.harris <- run.IV.county.Harris(MVTHEFT_countysumrate, county, 'items')
vehicle.county.value.harris <- run.IV.county.Harris(MVTHEFT_countysumrate, county, 'value')


# reduced form at county level
homicide.red.county.items.harris <- reduced.county.Harris(MURDER_countysumrate, county, 'items')
homicide.red.county.value.harris <- reduced.county.Harris(MURDER_countysumrate, county, 'value')
robbery.red.county.items.harris <- reduced.county.Harris(ROBBERY_countysumrate, county, 'items')
robbery.red.county.value.harris <- reduced.county.Harris(ROBBERY_countysumrate, county, 'value')
assault.red.county.items.harris <- reduced.county.Harris(AGASSLT_countysumrate, county, 'items')
assault.red.county.value.harris <- reduced.county.Harris(AGASSLT_countysumrate, county, 'value')
vehicle.red.county.items.harris <- reduced.county.Harris(MVTHEFT_countysumrate, county, 'items')
vehicle.red.county.value.harris <- reduced.county.Harris(MVTHEFT_countysumrate, county, 'value')


## Exact replication of Harris Results -------
homicide.items.harris <- run.IV.Harris(Amurders, harris, 'items')
homicide.value.harris <- run.IV.Harris(Amurders, harris, 'value')
robbery.items.harris <- run.IV.Harris(Arobbery, harris, 'items')
robbery.value.harris <- run.IV.Harris(Arobbery, harris, 'value')
assault.items.harris <- run.IV.Harris(Aassault, harris, 'items')
assault.value.harris <- run.IV.Harris(Aassault, harris, 'value')
vehicle.items.harris <- run.IV.Harris(Avehicle, harris, 'items')
vehicle.value.harris <- run.IV.Harris(Avehicle, harris, 'value')

homicide.red.items.harris <- reduced.Harris(Amurders, harris, 'items')
homicide.red.value.harris <- reduced.Harris(Amurders, harris, 'value')
robbery.red.items.harris <- reduced.Harris(Arobbery, harris, 'items')
robbery.red.value.harris <- reduced.Harris(Arobbery, harris, 'value')
assault.red.items.harris <- reduced.Harris(Aassault, harris, 'items')
assault.red.value.harris <- reduced.Harris(Aassault, harris, 'value')
vehicle.red.items.harris <- reduced.Harris(Avehicle, harris, 'items')
vehicle.red.value.harris <- reduced.Harris(Avehicle, harris, 'value')

#---------------------------------------------------------#
## First stage --------
#---------------------------------------------------------#
harris_agency_items_first <- summary(homicide.items.agency.Harris$stage1)$coefficients
harris_agency_value_first <- summary(homicide.value.agency.Harris$stage1)$coefficients
harris_county_items_first <- summary(homicide.county.items.harris$stage1)$coefficients
harris_county_value_first <- summary(homicide.county.value.harris$stage1)$coefficients
harris_items_first <- summary(homicide.items.harris$stage1)$coefficients
harris_value_first <- summary(homicide.value.harris$stage1)$coefficients
first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first,
                   harris_items_first,
                   harris_value_first)
HPBM_first_stage_results <- tibble(first_stage = map(first_list, tidy))

# Save all results ------
agency_items_list <- list(homicide.items.agency.Harris,
                          robbery.items.agency.Harris,
                          assault.items.agency.Harris,
                          vehicle.items.agency.Harris)
agency_value_list <- list(homicide.value.agency.Harris,
                          robbery.value.agency.Harris,
                          assault.value.agency.Harris,
                          vehicle.value.agency.Harris)

county_items_list <- list(homicide.county.items.harris,
                          robbery.county.items.harris,
                          assault.county.items.harris,
                          vehicle.county.items.harris)
county_value_list <- list(homicide.county.value.harris,
                          robbery.county.value.harris,
                          assault.county.value.harris,
                          vehicle.county.value.harris)


HPBM_agency_items_results <- tibble(map(agency_items_list, tidy))
HPBM_agency_value_results <- tibble(map(agency_value_list, tidy))
HPBM_county_items_results <- tibble(map(county_items_list, tidy))
HPBM_county_value_results <- tibble(map(county_value_list, tidy))

agency_red_items_list <- list(homicide.red.agency.items.harris,
                          robbery.red.agency.items.harris,
                          assault.red.agency.items.harris,
                          vehicle.red.agency.items.harris)
agency_red_value_list <- list(homicide.red.agency.value.harris,
                          robbery.red.agency.value.harris,
                          assault.red.agency.value.harris,
                          vehicle.red.agency.value.harris)

county_red_items_list <- list(homicide.red.county.items.harris,
                          robbery.red.county.items.harris,
                          assault.red.county.items.harris,
                          vehicle.red.county.items.harris)
county_red_value_list <- list(homicide.red.county.value.harris,
                          robbery.red.county.value.harris,
                          assault.red.county.value.harris,
                          vehicle.red.county.value.harris)


HPBM_agency_red_items_results <- tibble(map(agency_red_items_list, tidy))
HPBM_agency_red_value_results <- tibble(map(agency_red_value_list, tidy))
HPBM_county_red_items_results <- tibble(map(county_red_items_list, tidy))
HPBM_county_red_value_results <- tibble(map(county_red_value_list, tidy))


harris_items_list <- list(homicide.items.harris,
                          robbery.items.harris,
                          assault.items.harris,
                          vehicle.items.harris)
harris_value_list <- list(homicide.value.harris,
                          robbery.value.harris,
                          assault.value.harris,
                          vehicle.value.harris)

harris_red_items_list <- list(homicide.red.items.harris,
                          robbery.red.items.harris,
                          assault.red.items.harris,
                          vehicle.red.items.harris)
harris_red_value_list <- list(homicide.red.value.harris,
                          robbery.red.value.harris,
                          assault.red.value.harris,
                          vehicle.red.value.harris)

HPBM_items_results <- tibble(map(harris_items_list,tidy))
HPBM_value_results <- tibble(map(harris_value_list,tidy))
HPBM_red_items_results <- tibble(map(harris_red_items_list,tidy))
HPBM_red_value_results <- tibble(map(harris_red_value_list,tidy))

# save
if(controlled_subset==F & year_subset==F & city_outliers_eliminated==T & crime_problem_obs_dropped==F & scale_instruments==T){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM.RData')
}

# Save all results for texreg
HPBM_agency_items_results <- map(agency_items_list, texreg::extract)
HPBM_agency_value_results <- map(agency_value_list, texreg::extract)
HPBM_county_items_results <- map(county_items_list, texreg::extract)
HPBM_county_value_results <- map(county_value_list, texreg::extract)

HPBM_agency_red_items_results <- map(agency_red_items_list, texreg::extract)
HPBM_agency_red_value_results <- map(agency_red_value_list, texreg::extract)
HPBM_county_red_items_results <- map(county_red_items_list, texreg::extract)
HPBM_county_red_value_results <- map(county_red_value_list, texreg::extract)

HPBM_items_results <- map(harris_items_list, texreg::extract)
HPBM_value_results <- map(harris_value_list, texreg::extract)
HPBM_red_items_results <- map(harris_red_items_list, texreg::extract)
HPBM_red_value_results <- map(harris_red_value_list, texreg::extract)

harris_agency_items_first <- homicide.items.agency.Harris$stage1
harris_agency_value_first <- homicide.value.agency.Harris$stage1
harris_county_items_first <- homicide.county.items.harris$stage1
harris_county_value_first <- homicide.county.value.harris$stage1
harris_items_first <- homicide.items.harris$stage1
harris_value_first <- homicide.value.harris$stage1
first_list <- list(harris_agency_items_first,
                   harris_agency_value_first,
                   harris_county_items_first,
                   harris_county_value_first,
                   harris_items_first,
                   harris_value_first)
HPBM_first_stage_results <- map(first_list, texreg::extract)

# save
if(controlled_subset==F & year_subset==F & city_outliers_eliminated==T & crime_problem_obs_dropped==F & scale_instruments==T){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM-for-tables.RData')
}

if(controlled_subset==T){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM-for-tables-controlled.RData')
}

if(year_subset==T){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM-for-tables-year-subset.RData')
}

if(city_outliers_eliminated==F){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM-for-tables-with-outliers.RData')
}

if(crime_problem_obs_dropped==T){
  save(HPBM_agency_items_results,
       HPBM_agency_value_results,
       HPBM_county_items_results,
       HPBM_county_value_results,
       HPBM_agency_red_items_results,
       HPBM_agency_red_value_results,
       HPBM_county_red_items_results,
       HPBM_county_red_value_results,
       HPBM_items_results,
       HPBM_value_results,
       HPBM_red_items_results,
       HPBM_red_value_results,
       HPBM_first_stage_results,
       file = 'data/results-HPBM-for-tables-crime-obs-dropped.RData')
}



